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Abstract 

An  adaptive  grid  generation  method  is  presented 
which  is  based  on  the  reduction  of  truncation  error  in  the 
computational  plane.  Grid  adaption  is  achieved  through 
minimization  of  the  third  derivative  of  the  dependent 
variable  in  the  computational  plane.  Burger’s  equation  is 
solved  because  it  represents  typical  nonlinear  fluid-flow 
equations.  An  optimized  Successive-Over-Relaxation  method 
is  used  to  solve  Burger's  equation  using  second-order- 
upwind  differences  for  the  convective  term. 

Results  are  presented  for  several  cases  which  are 
compared  to  the  results  of  a  test  case.  The  results  show 
grid  points  are  concentrated  in  high  gradient  regions. 


Vlll 


ADAPTIVE  GRID  GENERATION  FOR  NUMERICAL 


SOLUTION  OF  BURGER'S  EQUATION 


I 


I. 


Introduction 


Finite-difference  algorithms  can  be  used  to  solve 
the  governing  equations  of  fluid  mechanics  and  heat  trans¬ 
fer  when  exact  solutions  are  unavailable.  In  order  to  get 
a  reasonable  solution  to  a  finite-difference  algorithm,  a 
suitable  grid  must  be  chosen.  In  most  algorithms,  a  grid 
is  selected  initially  and  the  grid  point  locations  are  held 
fixed  during  the  calculation.  The  best  grid  point  locations 
for  a  given  problem  are  determined  from  a  priori  knowledge 
of  the  exact  solution.  Numerical  grid  generation  schemes 
(see  1;  2)  can  be  used  to  determine  grid  point  locations 
without  a  priori  knowledge  of  the  exact  solution.  Adaptive 
grid  generation  schemes  (see  3;  4),  a  class  of  numerical 
grid  generation  schemes,  can  also  be  used  to  determine 
grid  point  locations;  however,  the  grid  point  locations 
are  not  held  fixed  during  the  calculation.  Movable  grid 
points  can  be  concentrated  in  high  gradient  regions,  a 
desirable  goal. 

The  purpose  of  this  thesis  is  to  develop  an  adap¬ 
tive  grid  method  which  reduces  the  truncation  error  in  the 
computational  plane.  The  method  will  be  based  on  the 


boundary  layer-dependent  grid  generation  method  of  Hodge 
and  Stone  (5) .  The  grid  spacing  is  determined  such  that 
the  truncation  error  is  minimized.  The  viscous  Burger's 
equation  will  be  used  as  a  test  of  the  adaptive  grid  method 


Background 

Numerical  grid  generation  schemes  (including  adap¬ 
tive  grid  generation  schemes)  should  satisfy  the  following: 

(1)  Grid  points  should  be  surface  oriented  so  boundary  con¬ 
ditions  can  be  applied  at  the  surface  without  requiring 
interpolation.  Thus,  interpolation  errors  are  not  intro¬ 
duced  into  the  solution.  In  an  adaptive  grid,  surface- 
oriented  grid  points  should  move  due  to  boundary  motion. 

(2)  High  gradient  regions,  such  as  shocks  and  boundary 
layers,  should  be  accurately  resolved  since  large  numerical 
errors  can  occur. 

A  surface-oriented  grid  can  be  generated  by  solving 
a  system  of  elliptic  partial  differential  equations,  one 
equation  for  each  coordinate  direction.  Thompson,  Thames, 
and  Mastin  (6)  develop  a  grid  generation  method  based  on 
the  solution  of  Poisson  equations.  Grid  point  placement  is 
controlled  by  the  forcing  functions  (also  called  the  grid 
control  functions)  in  the  Poisson  equations.  The  grid 
point  locations  are  found  by  interchanging  the  dependent 
and  independent  variables  in  the  Poisson  equations  and 
solving  the  resulting  equations. 


In  the  past,  gradients  were  resolved  by  placing 
many  grid  points  in  suspected  high  gradient  regions.  If 
the  locations  of  high  gradient  regions  were  unknown,  a  fine 
grid  was  used  over  the  entire  domain.  Since  more  points 
are  present  in  the  fine  grid,  more  computer  time  is  needed 
to  find  the  solution.  An  adaptive  grid  can  be  generated 
such  that  grid  points  will  move  into  high  gradient  regions 
as  the  gradients  in  the  computed  solution  develop.  Two 
benefits  accrue  from  the  use  of  an  adaptive  grid.  First, 
fewer  grid  points  are  needed  since  the  grid  points  will 
move  toward  high  gradient  regions,  thus  reducing  computer 
time.  Second,  truncation  error  in  the  computed  solution  is 
reduced.  Large  truncation  errors  can  occur  in  high  gradient 
regions.  If  grid  points  are  concentrated  in  high  gradient 
regions,  the  grid  spacing  is  reduced,  thereby  reducing 
truncation  error.  If  a  suitable  grid  adaption  criterion  is 
selected,  such  as  minimizing  the  magnitude  of  the  deriva¬ 
tives  in  the  truncation  error,  the  truncation  error  is 
reduced. 

Hodge  and  Stone  (5)  develop  a  boundary  layer- 
dependent  grid  based  on  the  Blasius  solution  which  mini¬ 
mizes  truncation  error.  Truncation  error  is  reduced  by 
placing  grid  points  such  that  increments  in  the  velocity 
are  constant.  Therefore,  the  second  and  higher-order 
derivatives  of  the  velocity  in  the  computational  plane 
which  contribute  to  truncation  error  are  minimized.  Grid 


point  placement  is  accomplished  by  modifying  the  grid  con¬ 
trol  functions  in  the  grid  generation  equations  of  Thompson 
et  al.  Since  the  Blasius  solution  was  used  to  place  the 
grid  points,  this  type  of  grid  generation  is  valid  only 
where  boundary  layer  theory  is  valid.  The  method  of  Hodge 
and  Stone  is  said  to  be  a  first-order  method.  The  grid 
generation  method  developed  in  this  thesis  is  similar 
except  for  the  following:  (1)  The  present  method  is  said  to 
be  second  order  because  the  third  derivative  of  the  velocity 
in  the  computational  plane  is  minimized.  (2)  The  grid  is 
not  dependent  on  boundary  layer  theory. 

Examples  of  adaptive  grid  methods  follow. 

Literature  Survey  of  Adaptive 
Grid  Methods 

Gradients  in  the  dependent  variables  can  be  used  to 
move  grid  points  as  the  computed  solution  develops.  Dwyer, 
Kee,  and  Sanders  (7)  place  grid  points  in  proportion  to 
gradients  that  appear  in  the  developing  computed  solution. 
Grid  skewness  is  one  weakness  of  this  method. 

Error  minimization  is  another  adaptive  grid  tech¬ 
nique.  Pierson  and  Kutler  (8)  minimize  a  measure  of  local 
truncation  error  with  respect  to  the  transformation  param¬ 
eters.  Anderson  and  Rai  (9)  induce  a  velocity  in  each  grid 
point  with  respect  to  every  other  grid  point  based  on  the 
gradient  at  each  grid  point  compared  to  the  average 
gradient.  Grid  points  are  attracted  to  regions  where  local 


gradients  are  higher  than  the  average  gradient  and  repelled 
from  regions  where  local  gradients  are  lower  than  the 
average.  Error  is  reduced  since  grid  points  are  concen¬ 
trated  in  high  gradient  regions.  Brown  (10)  minimizes 
truncation  error  in  the  computational  plane  by  generating 
a  grid  such  that  the  third  derivative  of  the  dependent 
variable  in  the  computational  plane  is  minimized. 

Variational  techniques  are  also  used  in  adaptive 
grid  generation.  Saltzman  and  Brackbill  (11)  utilize  a 
variational  technique  to  minimize  a  linear  combination  of 
grid  smoothness,  orthogonality,  and  area  variation. 

Other  criteria  can  be  used  to  adapt  a  grid  to  the 
developing  solution.  The  adaptive  grid  criterion  used  by 
Ghia,  Ghia,  and  Shim  (12)  is  the  minimization  of  the  coeffi 
cient  of  the  convective  term  in  the  transformed  flow  equa¬ 
tions.  This  technique  is  developed  for  high-Reynolds 
number  flows. 


Objectives 

The  objective  of  this  thesis  is  to  develop  an  adap¬ 
tive  grid  method  which  reduces  truncation  error  in  the 
computational  plane.  The  grid  control  function  used  in  the 
grid  generation  equation  is  systematically  determined  by 
minimizing  the  third  derivative  of  the  dependent  variable 
in  the  computational  plane.  The  method  is  based  on  the 
boundary  layer-dependent  method  of  Hodge  et  al.  and  is 
similar  to  the  method  of  Brown.  The  present  method  differs 


from  the  method  of  Brown  in  the  method  of  solving  the  grid 
generation  equation,  and  the  determination  of  the  grid  con¬ 
trol  function. 

A  one-dimensional  model  problem  (the  viscous 
Burger's  equation)  is  used  to  test  the  method.  Boundary 
conditions  are  imposed  such  that  a  boundary  layer-type 
velocity  profile  is  produced.  Grid  points  are  expected  to 
be  concentrated  in  the  boundary  layer. 

Overview 

Section  II  contains  the  following:  the  transforma¬ 
tion  used  to  transform  derivatives  in  the  physical  plane  to 
a  computational  plane;  the  one-dimensional  grid  generation 
equation;  the  criterion  used  to  adapt  the  grid;  the  deter¬ 
mination  of  the  grid  control  function;  the  one-dimensional 
model  problem;  the  solution  procedure;  and  the  error 
analysis  used  to  evaluate  the  method. 

The  criterion  used  to  adapt  the  grid  is  the  minimi¬ 
zation  of  the  third  derivative  of  the  velocity  in  the 
computational  plane.  Minimization  is  achieved  by  fitting 
the  computed  solution  with  a  second-order  polynomial.  The 
third  derivative  is  related  to  the  grid  control  function, 
thus  the  grid  points  are  moved  as  the  third  derivative  is 
minimized.  The  grid  control  function  is  determined  from  a 
Newton-Raphson  iterative  scheme  which  minimizes  the  third 
derivative . 
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Section  III  contains  the  discussion  of  results. 
Several  cases  were  run  to  discover  the  effect  of  various 
input  parameters  on  the  method.  In  most  cases,  grid  points 
were  concentrated  in  the  boundary  layer;  however,  trunca¬ 
tion  error  was  not  minimized  in  all  cases. 

The  conclusions  are  found  in  Section  IV.  The  method 
is  very  sensitive  to  initial  conditions;  therefore,  the 
method  may  only  be  useful  for  a  few  combinations  of  input 
parameters . 

Recommendations  for  further  work  include  using  a 
different  criterion  for  grid  convergence;  updating  the  grid 
periodically  instead  of  every  time  step;  and  using  a  differ¬ 
ent  method  for  determining  the  grid  control  function.  All 
recommendations  are  found  in  Section  V. 


II.  Mathematical  Formulation 


Transformation  from  Physical 
Plane  to  Computational  Plane 

Just  as  cylindrical  coordinates  are  used  instead  • 
of  cartesian  coordinates  to  simplify  the  solution  of  cer¬ 
tain  problems,  such  as  flow  around  a  cylinder,  transform¬ 
ing  variables  in  the  physical  plane  to  variables  in  a 
computational  plane  may  simplify  the  solution  of  a  numeri¬ 
cal  algorithm.  Arbitrarily-shaped  regions  in  the  physical 
plane  can  be  uniquely  transformed  to  a  rectangular  region 
in  the  computational  plane  with  uniform  grid  spacing, 
assuming  the  Jacobian  of  the  transformation  is  non-zero. 

The  shape  and  spacing  of  the  grid  in  the  computational  plane 
will  remain  constant  even  if  the  grid  points  are  moving  in 
the  physical  plane. 

For  a  one-dimensional  time-dependent  problem,  the 
computational  variables,  ?  and  t,  are  functionally  related 
to  the  physical  variables,  x  and  t,  by 

T  =  t  (la) 

C  =  C  (x,t)  (lb) 

Derivatives  in  the  physical  plane  can  be  transformed  in 
the  computational  plane  using  the  chain  rule  yielding 


where  the  subscripts  denote  partial  differentiation.  The 
dependent  flow  variable,  u,  can  represent  velocity,  tempera 
ture ,  pressure,  etc.  The  complexity  of  Eqs.  (2)  do  not 
outweigh  the  advantage  of  using  a  rectangular  grid  with 
uniform  spacing. 

Grid  Generation 

In  this  study,  a  one-dimensional  grid  is  used  for 
ease  of  computation  and  computer  economy. 

The  one-dimensional  grid  generation  equation  used 
in  this  thesis  is  developed  in  Brown  (10)  and  has  the  form 

+  Px^  =  0  (3) 

where  P,  the  grid  control  function,  is  a  function  of  r  and. 
t.  Brown  used  an  implicit  optimized  Successive-Over- 
Relaxation  (SOR)  method  to  solve  Eq.  (3).  A  second-order- 
upwind  difference  was  used  for  x  based  on  the  sign  of  P. 

A  central  difference  was  used  for  x  .  The  solution  of 
Eq.  (3)  using  the  given  finite  differences  produces 


2 

truncation  error  of  order  (PAc)  (13)  ,  where  Ac  is  the 
grid  spacing  in  the  computational  plane. 

A  more  accurate  method  of  solving  Eq.  (3)  is  given 
by  Lick  and  Gaskins  (13) .  A  Unified  Difference  Relation 
(UDR)  is  developed  which  is  an  exact  representation  of 
Eq.  (3)  for  a  constant  P.  The  UDR  corresponding  to  Eq.  (3) 
is  given  by 

xi+1  -  (l+e”P)xi  +  e“Pxi_1  =  0  (4a) 

If  P  takes  on  large  negative  values,  an  overflow  condition 
may  result  when  Eq.  (4a)  is  solved  on  a  computer.  Thus, 
for  P  negative,  the  UDR  corresponding  to  Eq.  (4a)  is  given 
by 

ePxi+1  -  (eP+l)xi  +  =  0  (4b) 

Equations  (4)  are  solved  using  a  tridiagonal  solver. 
Adaptive  Grid  Criterion 

Any  finite-difference  representation  of  a  deriva¬ 
tive  has  truncation  error  associated  with  it  due  to  the 
truncation  of  the  Taylor  series  used  to  define  the  deriva¬ 
tive.  For  example,  a  central-difference  representation  of 
a  first  derivative  in  the  computational  plane  is 

u  =  (ui+l"ui-l) /2Ac  +  T-E'  (5) 

where  the  truncation  error,  T.E.,  is  given  by 
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T.E.  =  - 


U44444  + 


Note  that  Eq.  (  5  )  is  said  to  be  second-order  accurate  and 
note  that  A  4  is  an  arbitrary  constant  and  for  convenience 
will  be  taken  as  unity  for  the  remainder  of  this  work. 

In  order  to  reduce  the  truncation  error  in  the 
computational  plane,  the  method  developed  in  this  thesis 
reduces  the  magnitude  of  u  . 

If  u  can  be  exactly  represented  by  a  second-degree 
polynomial  in  the  computational  plane,  u  and  all  higher- 

'-y 

order  derivatives  are  identically  zero.  Globally  approxi¬ 
mating  u  by  a  second-degree  polynomial  may  not  be  feasible. 
A  local  second-degree-polynomial  approximation  to  u  will 
minimize  u  locally,  thus  u  can  be  approximated  by 

S 


U(0  =  a2  4  +  aL  e  +  a 

1  11 


To  discover  what  values  of  u  contribute  to  the  local  trunca¬ 
tion  error,  approximate  u  by  a  central  difference,  thus 

4  4  Cf 


u  ~  u  -  u 

S«i  «i+1  Wi-i 


Using  central  differences  to  approximate  the  second  deriva¬ 
tives  in  Eq.  (8a)  yields 


Vq  '  Ui+2  -  2ui+l  +  2ui-l  -  ui-2 


Thus,  local  truncation  error  is,  in  a  general  way,  depen¬ 
dent  on  u  at  surrounding  grid  points.  This  suggests  that 
the  computational  plane  can  be  divided  up  into  several  sub- 
regions,  each  of  which  satisfies  Eq.  (7).  To  determine  the 

coefficients  a2  '  ai  '  anc^  ao  '  a  least-s<Juares  curve  fit 
i  1  ui 

is  performed  at  every  point  using  data  from  the  surrounding 
grid  points  at  the  previous  iteration  level.  The  least- 
squares  curve  fit  used  in  this  thesis  is  taken  from  an 
algorithm  in  Conte  and  de  Boor  (14).  The  optimum  curve  fit 
is  found  by  minimizing  a  sum  of  polynomials  of  degree  two 
and  below. 

Each  subregion  of  the  computational  plane  consists 

of  a  number  of  data  points,  N  ,  which  must  contain  at  least 

P 

three  points  and  not  more  than  the  total  number  of  grid 

points.  If  N  is  equal  to  the  total  number  of  grid  points, 
P 

only  one  set  of  coefficients  ,  ai  ,  and  ag  is  fQund 

l  i  i 

for  the  entire  computational  plane.  If  is  less  than  the 
total  number  of  grid  points,  a  balanced  curve  fit  is  per¬ 
formed  at  each  point;  i.e.,  if  N  =5,  the  values  of  u 

P 

used  to  calculate  the  curve  fit  will  be  taken  from  the 
i  ±  2,  i  ±  1,  and  i  locations.  At  or  near  the  boundaries, 
the  curve  fit  will  be  unbalanced. 

The  grid  control  function  must  now  be  determined 
such  that  the  grid  points  will  be  placed  in  such  a  way  as 
to  force  u  to  locally  fit  a  second-order  polynomial.  Thus 


u^  will  be  minimized  locally  and  truncation  error  will 
be  reduced. 

Grid  Control  Function  Evaluation 

Since  grid  point  locations  are  determined  from  the 

value  of  P  (see  Eqs.  (3)  and  (4)),  a  relation  is  developed 

involving  u  ,  P,  and  a„  which  will  be  minimized  in  a 
^^i  1 

Newton-Raphson  scheme.  Starting  with  an  equation  of  the 

form  of  Eq.  (8a),  u  can  be  approximated  by  a  difference 

i 

of  second  derivatives.  Thus 


u  ~  u  -  u 

1  X1  12 


(9) 


where  u  is  derived  from  Eq.  (2a)  and  u  is  derived 

X1  12 
from  Eq.  (7) .  From  Eq.  (2a) 


u  =  u 
?  x  C 


(10) 


Differentiating  Eq.  (10)  with  respect  to  5  yields 


u  =  u  x  +  u  x 
5C  xx  c  x  CC 


(11) 


Equation  (11)  is  used  to  define  u  .  To  determine  u  , 

X1  x2 

differentiate  Eq.  (7)  with  respect  to  q  twice,  thus 


UGC-  2a2  . 
x2  1 


(12) 
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Substituting  Eqs.  (11)  and  (12)  into  Eq.  (9)  yields  the 
equation  to  be  minimized: 


F  =  u  x^  +  u  “  2a_  =  0 
xx  £  x  ^  2 


To  get  F  as  a  function  of  P,  divide  Eq.  (13)  by  x  and 
note  x  /x  =  -P  from  Eq.  (3) ,  thus  F  becomes 


G(P)  =  u  x  -  Pu - -  =  0 

n  ^  x  x^ 


Values  of  P  can  be  found  by  using  a  Newton-Raphson  scheme 
which  has  the  form 


P(S+1)  =  p(S)  _  G(s)/(Gt)  (S) 


where  S  denotes  the  iteration  level  and  the  prime  denotes 


differentiation  with  respect  to  P.  Substituting  Eq.  (14) 
into  Eq.  (15)  gives 


p(S+U  =  p(S)  _  (u  x  _  Pu - —  )  /  (— u  )  (16) 

xx  c  x  x^  x 

Utilizing  Eqs.  (2)  for  u  and  u  and  after  some  simplifica- 

X  XX 


tion,  Eq.  (16)  becomes 


p(S+l)  =  p(S)  +  (u^  _  2a2)/u^  (17) 

Eq.  (17)  is  used  to  determine  P  at  each  grid  point.  Note 
if  P  is  converged,  the  second-degree-polynomial  approxima¬ 
tion  of  u  is  satisfied.  Thus,  the  optimum  grid  is  produced 
and  truncation  error  is  minimized. 


••  '•  A  •-  . 
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One-Dimen sional  Model  Problem 


To  test  the  adaptive  grid  method  given  above,  the 
one-dimensional  unsteady  viscous  Burger's  equation  is  used. 
The  viscous  Burger' s  equation  was  chosen  because  it  simu¬ 
lates  nonlinear  fluid  flow  equations.  Burger's  equation 
also  has  an  exact  solution  thus  the  truncation  error  can 
be  determined  directly. 

In  nonconservative  form,  the  viscous  Burger's  equa¬ 
tion  in  the  physical  plane  is 

ufc  +  uux  =  uxx/Re  (18a) 

where  Re  is  the  Reynolds  number.  The  boundary  conditions 
used  in  this  thesis  are 

u(-l,t)  =  1 

u (0  , t)  =  0  (18b) 

and  the  initial  conditions  are 

u  (0 , 0)  =  0 

u(x,0)  =1,  x  <  0  (18c) 

With  these  conditions,  a  velocity  profile  similar  to  a 
boundary  layer  profile  develops  near  the  x=0  boundary. 

The  thickness  of  this  boundary  layer-type  profile  is  pro¬ 
portional  to  1/Re.  The  exact  steady  state  solution  is 
given  by 

=  -  tanh(Re  ^) 


u  (x) 


(19) 


The  solid  line  in  Fig.  12  shows  u(x)  for  Re  =  1000.  Trans¬ 
forming  Eq.  (18a)  to  the  computational  plane  yields 


u  + 

T 


u-x 


u 


_L_  a, 

Re  Vc 


This  viscous  term  can  be  broken  into  two  parts, 
can  be  included  in  the  convective  term,  thus 


(20) 

one  of  which 


i 


I 


I 


u  + 

T 


u-x 


X 


JLL. 


Re  x 


3> 


u. 


JL£_ 


Re  x 


(21) 


Eq.  (21)  is  used  by  Brown  (10) .  Both  (20)  and  (21)  are 
used  in  this  thesis. 

An  implicit  optimized  SOR  method  is  utilized  to 
solve  Eqs.  (20)  and  (21) .  First-order-backward  differ¬ 
ences  are  used  for  and  the  grid  speed,  x^;  a  second- 
order-upwind  difference  is  used  for  u  based  on  the  sign 
of  its  coefficient;  central  differences  are  used  for  the 
viscous  terms  and  x^ . 


Solution  Procedure 

A  listing  of  the  computer  program  BURG  used  to 
solve  the  viscous  Burger's  equation  using  an  adaptive  grid 
is  found  in  Appendix  B.  The  algorithm  used  in  the  program 
consists  of  the  following  steps: 

Step  1.  Input  an  initial  guess  for  the  grid  con¬ 
trol  function,  P. 
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Step  2.  Solve  the  appropriate  grid  generation 
equation,  either  Eq.  {3)  or  Eq.  (4)  for  the  grid  point  loca¬ 


tions  . 

Step  3.  Solve  the  appropriate  equation,  either 
Eq.  (20)  or  Eq.  (21),  for  the  velocity,  u,  at  each  grid 
point . 

Step  4.  Perform  a  least-squares  curve  fit  of  the 

computed  solution  at  every  grid  point  to  determine  a 2  . 

i 

Step  5.  Using  Eq.  (17)  ,  generate  the  values  of  the 
grid  control  function,  P,  at  every  grid  point. 

Step  6.  Repeat  steps  2  through  5  until  steady 
state  is  reached.  Steady  state  is  assumed  to  be  reached 
when  the  average  difference  in  the  computed  solution  between 
two  time  levels  is  less  than  a  specified  tolerance;  i.e.. 


1 

I 


I 

I 

i=l 


n 
u . 
i 


u 


n-1 

i 


<  £ 


(22) 


where  I  is  the  total  number  of  grid  points,  n  denotes  the 
time  level,  and  the  tolerance,  e,  has  a  value  of  10  ®  in 
this  thesis. 


Error  Analysis 

Three  types  of  error  can  be  used  to  evaluate  the 
success  of  this  method:  truncation  error,  fit  error,  and 
steady  state  error.  Since  the  exact  solution  is  known. 


truncation  error  can  be  calculated  as  the  difference  between 


the  computed  solution  and  the  exact  solution.  Fit  error  is 
defined  as  the  difference  between  the  computed  solution  and 
the  least-squares  curve  fit.  Average  values  of  truncation 
error,  T.E.  ,  and  fit  error,  F.E.  ,  are  calculated  as 

dV  6  a.  v  G 


T.E. 


ave 


.1  l 

1  i-i 


u .  - 

l 


u 


exact 1 


(23a) 


!  I 

F.E.  =  ±  Z  I u.  - 
ave  I  .  . 1  l 


u 


fit1 


(23b) 


Steady  state  error  is  defined  as  the  difference  in  the 
computed  solution  between  the  final  two  time  levels;  i.e.. 


S.S.E. 


(23c) 


where  nT  denotes  the  time  level  at  convergence.  Steady 

state  error  gives  an  indication  of  local  convergence. 

The  optimum  grid  is  defined  when  the  average  grid 

speed  is  less  than  a  specified  tolerance,  x  ,  thus 

Tacc 


i  I 

T  ^  |x  |  <  x  (23d) 

i=l  Ti  Tacc 


When  condition  (23d)  is  met,  the  grid  is  held  fixed  and  the 
computed  solution  is  allowed  to  iterate  until  condition 
(22)  is  met.  When  the  computed  solution  has  converged,  the 
fit  error  is  assumed  minimized,  thus  the  truncation  error 


is  assumed  minimized. 


III.  Discussion  of  Results 


To  test  the  adaptive  grid  method  developed  in  this 
thesis,  several  cases  were  run  on  a  CD C  Cyber  845  computer 
to  determine  the  effects  of  various  input  parameters. 

The  conditions  used  for  the  different  cases  are  given  in 
Table  1.  Results  are  given  in  Table  2  and  also  in  Figures 
1  through  32.  (All  figures  are  found  in  Appendix  A.) 
Figures  1  through  10  show  the  initial  and  final  grids.  An 
idea  of  how  well  the  grid  adapted  can  be  seen  in  these 
figures.  Figures  11  through  21  show  the  exact  and  computed 
velocity  profiles.  How  well  the  method  matched  the  exact 
solution  can  be  seen  in  these  figures.  Figures  22  through 
32  show  the  steady  state  error  at  each  grid  point  in  the 
computational  plane. 

Test  Case  (Case  1J 

The  results  of  a  test  case  (Case  1)  are  used  for 

comparisons  with  the  results  of  the  other  cases.  The 

parameters  are  Reynolds  number.  Re;  average  grid  speed 

tolerance,  xT  ;  minimum  u  derivative  allowed,  u^;  the 
acc 

number  of  points  used  in  the  least- squares  curve  fit,  N  ; 

P 

the  number  of  grid  points  in  which  u  is  approximated  by  a 
linear  function.  Also,  the  values  of  the  grid  control 
function  P;  the  total  number  of  grid  points,  I;  and  the 
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CONDITIONS  FOR  CASES  1  THROUGH  11 


11  does  not  use  an  adaptive  grid. 


TABLE  2 


RESULTS  FOR 

CASES  1 

THROUGH  11 

Case 

T.E. 

ave 

F.E. 

ave 

Total  No. 

of  Time 
Steps, 

No.  of  Points 
in  the  Boundary 
Layer 

Total 
CPU  Time 
(seconds) 

1 

.01168 

.05590 

8 

6 

1.082 

2 

.002985 

.01935 

26 

11 

4.930 

3 

.04158 

.1336 

22 

7 

4.074 

4 

.001389 

. 003228 

45 

8 

8.328 

5 

.001024 

.0005713 

7 

9 

3.166 

6 

.004002 

.02276 

18 

6 

2.400 

7 

.004687 

.01895 

10 

6 

1.500 

8 

. 004127 

.06411 

11 

7 

1.665 

9 

.009522 

.05132 

7 

5 

1.158 

10 

.001911 

.02685 

7 

2 

.823 

11 


002021 


3 


2 
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equations  used  to  generate  the  grid  and  the  form  of  the 

viscous  Burger's  equation  are  investigated. 

Reynolds  number,  Re,  is  varied  to  determine  how  well 

the  grid  adapts  as  the  boundary  layer  thickness  changes  and 

the  test-case  value  of  Re  is  1000. 

The  value  of  the  grid  convergence  criterion,  x  , 

Tacc 

is  varied  to  determine  whether  the  optimum  grid  is  found. 

Fit  error  is  expected  to  be  minimized  for  a  small  value  of 
x  and  the  test-case  value  is  .005. 

T 

acc 


The  value  of  u^  must  be  limited  to  non- zero  values 

since  u^  is  in  the  denominator  of  the  second  term  on  the 

right-hand  side  of  Eq.  (17) .  Test-case  value  of  u^  is 

chosen  as  .1  because  u^  =  .1  represents  ten  points  in  the 

boundary  layer  if  u  is  assumed  linear  there. 

The  number  of  points  used  in  the  least- squares 

curve  fit,  N  ,  is  varied  to  determine  the  dependence  of 
P 

truncation  error  on  the  computed  solution  at  the  surround¬ 
ing  grid  points.  The  test-case  value  of  was  chosen  as 
5  based  on  Eq.  (8b) . 

Early  results  showed  that  very  few  grid  points  were 

concentrated  in  the  boundary  layer.  Positive  values  of  the 

grid  control  function,  P,  are  needed  to  move  points  into 

the  boundary  layer.  In  order  to  force  P  to  be  positive  in 

the  boundary  layer,  u  is  approximated  as  a  linear  function 

in  a  region  next  to  the  x=0  boundary.  In  this  region,  a 2 

c  i 

is  zero,  thus  P  will  become  positive  if  u  is  negative. 
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since  u^  is  negative  in  the  boundary  layer  (see  Eq.  (17) ) . 
Note  that  u  is  still  minimized  if  u  is  approximated  by 
a  linear  function.  The  number  of  points  in  the  linear 
r-  gion  is  five  in  the  test  case. 

If  P  changes  sign  somewhere  in  the  domain,  a  rela¬ 
tively  large  spacing  is  produced  between  two  grid  points. 
This  spacing  is  a  result  of  a  negative  P  causing  grid 
points  to  move  away  from  the  boundary  layer  and  a  positive 
P  causes  grid  points  to  move  toward  the  boundary  layer. 

This  spacing  is  undesirable  because  a  smooth  transition  is 
needed  between  regions  of  high  gradient  and  regions  of  low 
gradient  for  an  accurate  solution  (10).  To  eliminate  the 
large  spacing,  P  can  be  constrained  to  non-negative  values. 
In  the  test  case,  P  is  allowed  to  take  on  negative  values. 

The  total  number  of  grid  points,  I,  is  varied  to 
determine  if  grid  refinement  leads  to  a  more  accurate  solu¬ 
tion.  A  larger  number  of  grid  points  in  the  boundary  layer 
is  expected  as  the  grid  is  refined  and  the  test-case  value 
of  I  is  21. 

The  form  of  the  equations  used  or  the  method  used 
to  solve  them  may  have  an  effect  on  the  solution.  Equa¬ 
tions  (4)  are  used  for  grid  generation  in  the  test  case,  and 
Eq.  (20)  is  used  to  solve  the  viscous  Burger's  equation. 

As  a  comparison,  an  optimized  SOR  method  is  used  to  solve 
both  Eq.  (3),  for  the  grid  point  locations,  and  Eq.  (21), 
for  the  solution  of  Burger's  equation. 


An  initial  guess  for  P  at  every  point  is  needed  to 
start  the  program.  An  initial  P  of  zero  produces  a  car¬ 
tesian  grid.  However,  early  results  showed  that  the  grid 
did  adapt  very  well  using  an  initial  P  of  zero.  Thus,  a 
positive  value  of  P  was  used  initially  to  produce  a  grid 
which  tends  toward  the  boundary  layer. 

An  error  comparison  is  made  between  the  test  case 
and  a  fixed-grid  case.  A  fixed  grid  is  achieved  by  holding 
P  constant  at  each  grid  point.  The  fixed  grid  corresponds 
to  the  initial  grid  of  the  test  case. 

Expected  Results 

The  computed  solution  is  expected  to  match  the  exact 
solution,  thus  truncation  error  is  reduced,  and  grid  points 
are  expected  to  concentrate  in  the  boundary  layer.  As  a 
rule  of  thumb,  five  to  ten  grid  points  are  necessary  to 
accurately  resolve  the  gradient. 

Results  and  discussion  of  the  cases  tested  follow. 

Grid  Dependence  on  Re 

Case  8  was  used  to  determine  whether  more  grid 
points  can  be  forced  into  the  boundary  layer  if  a  small  Re 
is  used.  The  Reynolds  number  for  Case  8  was  100  and  seven 
points  were  placed  in  the  boundary  layer  as  compared  to 
six  for  the  test  case  (Fig.  19) .  Average  truncation  error 
was  less  than  half  the  average  truncation  error  of  the 
test  case  (see  Table  2) ,  due  to  a  smaller  gradient.  The 


Reynolds  number  for  Case  9  was  10,000  and  this  Reynolds 
number  was  used  to  see  the  effect  of  using  a  larger  gradient 
on  the  grid  adaption.  Five  points  were  placed  in  the  boun¬ 
dary  layer  (Fig.  20)  and  truncation  error  was  not  signifi¬ 
cantly  reduced.  As  Reynolds  number  is  increased,  more 
grid  points  are  expected  to  be  concentrated  in  the  boundary 
layer;  however,  for  the  grid  adaption  used,  fewer  grid 
points  were  concentrated  in  the  boundary  layer  as  Re  is 
increased. 


Grid  Dependence  on  x 


acc 

The  value  of  x 


in  Case  4  was  .0001  and  more 


time  steps  were  required  for  a  converged  solution  than  for 
the  test  case.  Average  fit  error  was  less  than  one- tenth 
the  fit  error  for  Case  1.  This  fact  suggests  a  more 
optimum  grid  was  produced  in  Case  4.  However,  the  grid 
produced  contained  a  very  large  grid  spacing  between  two 
points  (Fig.  4) .  The  velocity  profile  (Fig.  15)  shows  that 
the  curve  fit  is  not  very  good  in  the  vicinity  of  the  large 
grid  spacing.  Thus,  the  third  derivative  is  not  minimized 
near  the  large  spacing  and  truncation  error  is  not  mini¬ 
mized  there. 

The  assumption  that  a  converged  grid  leads  to 


minimum  truncation  error  is  not  correct. 


Grid  Dependence  on 

In  an  attempt  to  allow  P  to  be  determined  by  actual 
slopes  rather  than  limited  slopes,  the  value  of  u^  in  Case  6 
was  reduced  to  .01.  Very  few  values  of  u^  are  above  .1 
so  reducing  the  limit  on  u^  allows  more  values  of  the  slope 
to  affect  the  solution.  The  final  grid  produced  for  Case  6 
(Fig.  6)  had  an  extremely  large  grid  spacing  between  two 
points  which  caused  the  computed  solution  to  converge  very 
slowly.  Note  that  the  curve  fit  is  not  very  good  near  the 
large  spacing,  thus  truncation  error  is  not  minimized  there 
(Fig.  17). 

Grid  Dependence  on 

Poor  results  were  achieved  when  N  was  allowed  to 

P 

include  all  grid  points  (Case  3) .  One  of  the  largest  values 
of  average  truncation  error  was  produced  by  Case  3.  The 
results  show  a  parabola  is  not  a  good  fit  to  a  hyperbolic 
tangent. 

Grid  Dependence  on  Linear  Region 

The  linear  region  had  an  effect  on  the  shape  of 
the  final  grid  in  all  cases  where  a  linear  approximation 
for  u  is  made  (Figs.  1,  2,  4,  6,  7,  8,  10).  In  all  these 
cases,  the  final  grids  exhibit  a  relatively  large  grid 
spacing  between  two  grid  points.  These  spacings  occur  at 
or  near  the  edge  of  the  linear  region.  This  suggests  that 
may  be  causing  P  to  change  signs  since  a ^  is  zero  in 
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the  linear  region  and  non-zero  outside  the  region.  Because 
of  the  large  spacing,  a  good  curve  fit  cannot  be  performed, 
thus  truncation  error  is  not  minimized  (Figs.  12,  13,  15, 

17,  18,  20,  21). 

In  Case  2,  u  was  approximated  as  a  linear  function 
in  a  10-point  region  near  the  x=0  boundary.  The  number  of 
points  in  the  boundary  layer  for  Case  2  is  11.  Increasing 
the  size  of  the  linear  region  is  a  good  method  of  increasing 
the  number  of  points  in  the  boundary  layer;  however,  a  good 
curve  fit  is  not  produced  at  the  edge  of  the  linear  region 
(Fig.  2). 

Grid  Dependence  on  Non-Negative  P 

As  noted  above,  a  large  spacing  in  the  grid  is  pro¬ 
duced  near  the  edge  of  the  linear  region  where  P  changes 
sign.  To  eliminate  this  large  spacing,  P  was  limited  to 
non-negative  values  in  Case  7.  A  better  grid  is  produced 
(Fig.  7)  but  the  curve  fit  is  not  very  good  at  the  edge  of 
the  linear  region  (Fig.  18)  thus  truncation  error  is  not 
minimized. 

Effect  of  Grid  Refinement 

When  the  total  number  of  grid  points,  I,  is  allowed 
to  increase,  better  resolution  of  the  gradient  is  expected; 
i.e.,  a  proportionate  number  of  grid  points  is  expected  to 
move  into  the  boundary  layer.  In  Case  5,  I  is  51  and  the 
results  show  lower  truncation  error  and  fit  error,  but  only 
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nine  grid  points  were  placed  in  the  boundary  layer.  Twelve 
to  15  grid  points  were  expected  to  move  into  the  boundary 
layer  since  more  than  twice  the  number  of  grid  points  were 
available  for  adaption.  Not  much  adaption  took  place  in 
Case  5  (Fig.  5) ,  thus  similar  results  could  be  expected 
from  a  fixed  grid  with  the  same  number  of  points.  Note, 
however,  the  increase  in  CPU  time  over  Case  1. 

Effect  of  Equations  Used 

In  Case  10,  Eq.  (3)  was  used  to  generate  the  grid 
and  Eq.  (21)  was  used  to  solve  the  viscous  Burger's  equa¬ 
tion.  Both  were  solved  implicitly  using  an  optimized  SOE 
method  and  poor  results  were  found.  Only  two  points  were 
placed  in  the  boundary  layer.  Although  average  truncation 
error  is  low,  truncation  error  is  not  minimized  for  the 
points  in  the  boundary  layer  (Fig.  21) . 

Effect  of  Fixed  Grid 

A  constant  value  of  P  was  used  in  Case  11  to  hold 
grid  points  fixed.  The  results  were  poor  when  compared  to 
adaptive  grids.  The  average  truncation  error  was  low 
because  only  two  points  were  in  the  boundary  layer  where 
the  largest  truncation  errors  occur. 

Steady  State  Error 

Steady  state  error  is  defined  by  Eq.  (24d) .  Figures 
22  through  31  show  steady  error  in  the  computational  plane. 
The  regions  of  highest  error  correspond  to  the  linear 
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regions  in  most  cases.  This  suggests  that  the  computed 
solution  is  not  fully  converged  in  the  linear  region.  A 
nearly  uniform  distribution  of  error  was  expected. 

Figure  32  shows  a  comparison  of  steady  state  error 
between  Cases  11  and  1.  Note  that  for  the  constant  grid 
case  (Case  11) ,  only  two  points  were  in  the  boundary  layer 
while  six  points  were  in  the  boundary  layer  in  the  adaptive 
case  (Case  1).  Thus,  even  though  the  fixed  grid  case  has 
less  error,  the  boundary  layer  is  not  resolved. 


IV.  Conclusions 


The  following  conclusions  can  be  drawn  from  the 
results  presented: 

1.  The  adaptive  grid  method  concentrates  grid 
points  in  the  boundary  layer.  For  Reynolds  number  1000, 
an  average  of  about  seven  grid  points  is  placed  in  the 
boundary  layer. 

2.  Truncation  error  is  not  minimized  everywhere  in 
the  computational  plane  because  an  optimum  grid  is  not 
produced. 

3.  The  method  is  very  sensitive  to  initial  condi¬ 
tions.  A  change  in  one  input  parameter  produces  drastic 
changes  in  grid  shape  and  velocity  profile. 

4.  Final  grid  spacing  is  affected  by  the  linear 
region  where  u  is  approximated  by  a  linear  function.  A 
relatively  large  spacing  between  two  grid  points  is  pro¬ 
duced  at  the  edge  of  the  linear  region  because  the  grid 
control  function,  P,  changes  sign  there. 


V.  Recommendations 


Several  changes  to  the  method  are  recommended  in 
order  to  produce  an  optimum  grid  which  minimizes  truncation 
error. 

The  grid  is  assumed  converged  when  the  average  grid 

speed  is  less  than  x  .  A  better  convergence  criterion 

Tacc 

is  the  fit  error.  Determination  of  the  grid  control  func¬ 
tion,  P,  is  based  on  the  least- squares  curve  fit,  thus  if 
P  is  converged,  then  u  satisfies  the  second-degree-polynomial 
approximation  and  fit  error  is  minimized.  A  truly  optimum 
grid  is  produced  when  P  is  converged. 

In  the  present  method,  the  grid  is  updated  every 
time  step.  Better  adaption  may  be  achieved  by  periodically 
updating  the  grid.  Updating  in  the  grid  is  similar  to 
finding  the  computed  solution  on  a  fixed  grid,  then  using 
that  solution  as  an  initial  guess  for  the  next  solution. 

The  linear  region  affects  the  final  grid  spacing. 
However,  approximating  u  as  a  linear  function  near  the  x=0 
boundary  is  necessary  to  attract  grid  points  to  the  boun¬ 
dary  layer.  When  sufficient  grid  points  are  in  the  boundary 
layer,  a  linear  approximation  is  no  longer  needed,  thus  the 
grid  can  be  allowed  to  smoothe  itself  and  allow  better  curve 
fits  to  be  performed. 


Average  error  is  used  to  determine  if  the  computed 
solution  is  converged  (see  condition  (22) ) .  Using  maximum 
error  as  a  convergence  criterion  may  give  better  results 
and  a  more  uniform  distribution  of  steady  state  error. 

Lastly,  another  method  to  determine  the  grid  control 
function,  P,  may  be  necessary.  The  Newton- Raphson  used  in 
this  thesis  can  only  minimize  the  third  derivative  of  u  in 
the  computational  plane.  A  method  which  eliminates  the 
third  derivative  may  reduce  truncation  error  significantly. 


Appendix  A:  Figures 
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Figure  1 
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Figure  4 
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INITIAL  ANO  FINAL  GRID  FOR  CASE  9  (RE=10000) 

Figure  9 


Figure  11 
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VELOCITY  PROFILE  FOR  CRSE  2  (  10  PT  LINEAR  REGION) 
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Figure  13 


VELOCITY  PROFILE  FOR  CASE  3  INP=21J 

Figure  14 
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ERROR  BETWEEN  TINE  STEPS  FOR  CASE  1 

Figure  22 
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Figure  26 
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ERROR  BETWEEN  TIME  STEPS  FOR  CRSE  10 


Figure  31 
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